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Abstract 


Bursting neurons fire rapid sequences of action potential spikes followed by a quiescent period. 
The basic dynamical mechanism of bursting is the slow currents that modulate a fast spiking 
activity caused by rapid ionic currents. Minimal models of bursting neurons must include both 
effects. We considered one of these models and its relation with a generalized Kuramoto model, 
thanks to the definition of a geometrical phase for bursting and a corresponding frequency. We 
considered neuronal networks with different connection topologies and investigated the transition 
from a non-synchronized to a partially phase-synchronized state as the coupling strength is varied. 
The numerically determined critical coupling strength value for this transition to occur is compared 
with theoretical results valid for the generalized Kuramoto model. 
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I. INTRODUCTION 


Neurons are known to exhibit a plethora of dynamical behaviors, represented by the 
generation of action potential patterns. One of such patterns is bursting, defined by the 
repeated firing of action potentials followed by quiescent periods. Hence the dynamics of 
bursting neurons has two timescales: a fast scale related to spiking and a slow scale of 
bursting itself. These timescales are related to different biophysical mechanisms occuring at 
the level of neuron membrane: there are fast ionic currents (chiefly Na + and K + ) responsible 
for spiking activity and slower Ca ++ currents that modulate this activity. 

Most neurons exhibit bursting behavior if conveniently stimulated. For example, in the 
neocortical layer 5 pyramidal neurons, when stimulated with DC current pulses, fire an initial 
burst of spikes followed by shorter bursts [l, 3] . In layers 2, 3, and 4 chattering neurons 
fire high-frequency bursts of 3 — 5 spikes with a short interburst period .3,4]. Cortical 
interneurons have been found to exhibit bursting as a response to DC pulses 5 . Pyramidal 
neurons in the CAl region of hippocampus produce high-frequency bursts after current 
injection [(3]. Thalamocortical neurons and reticular thalamic nucleus inhibitory neurons 


exhibit bursting as weli [7], Pnrkinje ceils in cerebe.iunr can burst when their synaptic input 
is blocked [8[. Bursting is also an important feature of sensory systems, because bursts can 
increase the reliability of synaptic transmission { 9 ]. I11 some systems, bursts improve the 
signal-to-noise ratio_of sensory responses and might be involved in the detection of specific 
stimulus features (iq| . 

Due to both synaptic coupling and common inputs among neurons there are many types 
of synchronization, which can be generally regarded as the presence of a consistent temporal 


relationship between their activity patterns [11-13]. A strong form of the latter relationship 


is complete synchronization, where neurons spikes at the same time, i.e. a precise temporal 
coincidence of events. A weaker relationship is bursting synchronization, in which only the 
beginning of bursting is required to occur at the same time, even though the repeated spiking 
may not occur synchronously. 

There has been observed bursting synchronization in cell cultures of cortical neurons, 
where uncorrelated firing appeared withing the first three days and transformed progressively 


into synchronized bursting within a week 


15]. Large-scale bursting synchonization in the 


7 — lAHz range has been found in the thalamus during slow-wave sleep, partially originated 
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in the thalamus and gated by modulatory input from the brainstem 16]. Various areas of 
the basal ganglia have been found to exhibit bursting synchronization related to Parkinson’s 
disease and resting tremor 17]. 

There exists sound neurophysiological evidence that hypokinetic motor symptoms of 
Parkinsons disease such as slowness and rigidity of voluntary movements are closely related 
to synchronized bursting in the 10 — 30 Hz range [l8 21]. The connection between bursting 
synchronization and pathological conditions like Parkinson’s disease, essential tremor and 
epilepsy has led to the proposal of many control strategies aiming to suppress or mitigate 


bursting synchronization 


22 ], 


One of such strategies is deep-brain stimulation (DBS), which consists of the application 
of an external high-frequency (> 100 Hz) electrical signal by depth electrodes implanted in 
target areas of the brain like the thalamic ventralis intermedius nucleus or the subthalamic 


nucleus 


23]. The effect of DBS would be similar to that produced by tissue lesioning and 


has proved to be effective in suppression of the activity of the pacemaker-like cluster of 
synchronously firing neurons, and achieving a suppression of the peripheral tremor 24], 
There is strong clinical evidence that DBS is a highly effective technique for treatment of 

261 . 


patients with Parkinson’s disease 25 


In spite of these results, DBS is yet far from being completely understood. Many re¬ 
sults in this field have been obtained from empirical observations made during stereotaxic 
neurosur gery , but further progress can be obtained with proper mathematical modelling of 


DBS 


27 


29]. The effects of DBS in networks of bursting neurons have been investigated 


when DBS is implemented through an harmonic external current 30] and a delayed feedback 


signal 


ail- 


On modelling the response of a neuronal network to an external perturbation like DBS it is 
of paramount importance to keep the model simple enough such that large-scale simulations 
(using a large number of neurons) can be performed in a reasonable computer time. In such 
reductionist point of view a minimal model could be one in which we can assign a geometrical 
phase to the bursting activity. The bursting neuron is thus regarded as a phase oscillator 
undergoing spontaneous oscillations with a given frequency. Thus bursting synchronization 
becomes a special case of phase synchronization, a phenomenon well understood for coupled 


oscillators with and without external excitation 


32]. 


A simple model for the dynamics of nonlinearly coupled phase oscillators is the Kuramoto 


3 

















model, which in its original version considers a global (all-to-all) coupling 33]. It can be 
generalized by considering an arbitrary coupling architecture (generalized Kuramoto model) 


34], The particular interest in such models is that many analytical and numerical results are 


known for them, specially the global case for which a mean-held theory exists for the tran¬ 
sition between a non-synchronized to a (phase-)synchronized behavior 35]. For generalized 
Kuramoto models it is possible to derive analytical expressions for the critical value of the 


coupling strength for which the abovementioned transition occurs 


Hence such a body 


of knowledge can be applied to networks of bursting neurons, helping to design strategies of 
synchronization control and/or suppression like DBS. 


The main goal of this paper is to show, using analytical and numerical arguments, that 
a system of coupled bursting neurons described by Rulkov’s model can be reduced to a 
generalized Kuramoto model. This reduction is valid as long as phase synchronization is 
concerned, since for frequency synchronization the behaviors can be quite different, though. 
We consider, in particular, some widely used connection topologies, like random (Erdos- 
Renyi), small-world, and scale-free networks. We show that the analytical results for the 
critical coupling strength to synchronized behavior, originally derived for the generalized 
Kuramoto model, can be used to describe the synchronization transition also for networks 
of bursting neurons. 


As a matter of fact, since bursting activity presents two timescales it can be also ap¬ 
proached from the point of view of a relaxation oscillator j^]. In our work, however, we 
describe bursting using a single phase. This simplification is justified since phase synchro¬ 
nization of bursting is chiefly related to the slow timescale. In other words, the fast spikes 
can be nonsynchronized even though the slow dynamics is synchronized. 


This paper is organized as follows: in Section II we describe the model we used to describe 
bursting neurons. Section III considers networks of coupled bursting neurons using different 
connection topologies. Section IV reviews some results on the generalized Kuramoto model, 
and Section V includes the comparisons we made between Kuramoto model and the network 
of bursting neurons. Our Conclusions are left to the final Section. 
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II. MODELS OF BURSTING NEURONS 


The choose of a suitable model describing the dynamics of biological neurons is dictated 
by some requirements. First the model must take into account the kind of dynamics one 

For example, if all one needs is to describe a spiking neuron, 


wishes to describe 


for the sake of neural coding simulations for example, a simple leaky integrate-and-fire 
(LIF) model would be enough 39] • However if one needs to describe the interplay between 
different ionic currents flowing through the neuron membrane, the Hodgkin-Huxley (HH) 


model would be a natural choice 


40]. On the other hand, the HH model would require far 


more computational power than the LIF since the former involves four complicated first- 
order differential equations whereas the latter just one simple equation. 

With bursting neurons this criterion also holds. Given that bursting results from the 
interplay between fast and slow ionic currents, Hodking-Huxley-type models would need at 


least one more equation to describe slow Ca modulation 


41 


42], A mode 


sensitive neurons exhibiting bursting has been proposed by Huber and Braun 


43 


thermally 


45], which 


describes spike train patterns experimentally observed in facial cold receptors and hypotha¬ 


lamic neurons of the rat 


46], electro-receptors organs of freshwater catfish |47j], and caudal 


photo-receptor of the crayfish 


. However, the Huber-Braun model has 5 differential equa¬ 
tions for each neuron, and computational limitations impose restrictions to its use for large 
networks 0 . 

If numerical simulations do not need to take into account the effect of system parameters 
and only the phenomenological aspects of bursting are relevant, then a good choice is the 
two-dimensional mapping equations proposed by Rulkov 


x(n + 1) = 


a 


+ y{n), 


1 + [ x(n )] 2 
y(n + 1) = y(n ) — ax{n) — /3, 


( 1 ) 

( 2 ) 


where x is the fast and y is the slow dynamical variable, whose values are taken at discrete 
time t = nr, with r = 1 and n — 0,1,2 ,.... 

The parameter a affects directly the spiking timescale, its values being chosen in the 
interval [4.1,4.3] so as to produce chaotic behavior for the evolution of the fast variable x n , 
characterized by an irregular sequence of spikes. The parameters a and /?, on their hand, 
describe the slow timescale represented by the bursts, and take on small values (namely 
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FIG. 1. Time evolution of the (a) fast and (b) slow variables in the Rulkov map CO)-© for a = 4.1, 
a = /3 = 0.001. 


0 .001) so as to model the action of an external dc bias current and the synaptic inputs on 
a given isolated neuron 53]. 

The Rulkov model was derived using dynamical rather than biophysical hypotheses. We 
choose the parameter a so as to yield chaotic behavior for the characteristic spiking of the 
fast variable x(ri) [Fig. Ufa)]. The bursting timescale, on the other hand, comes about the 
influence of the slow variable y(n), which provide a modulation of the spiking activity due 
to a saddle-node bifurcation [Fig. CEfb)]. On comparing the dynamics of the Rulkov map 
with similar results of the Huber-Braun model one is led to an approximate correspondence 
between variables of the Rulkov map and variables with biophysical significance: the discrete 
time n in the map corresponds to 0.5ms of the continuous time; and x stands for the 
membrane potential, in such a way that each unit of x in Fig. Ufa) corresponds to circa 


20 mV. The Ru 
models 


49 


54 


kov model has been used in several numerical investigation of coupled neuron 


56] 


Due to the chaoticity of the ^-dynamics of the Rulkov map, the duration of each burst 
suffers a slight variability. The beginning of each burst can be chosen rather arbitrarily, but 
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it turns out that a useful choice is to consider the local maxima of the variable y: let ni 
be the time at which the £th burst begins [see Fig. [0(b)]. Hence the duration of this burst 
is ri£ + i — 7i£. We can define a geometric bursting phase by considering a variable ip that 
increases of 2ir after a bursting event. A linear interpolation gives 


p>(n) = 27: £ + 2tt 


n — ri£ 
n e+i — n t 


(' ni<n < ri£ +1 ). 


( 3 ) 


We can also define a bursting (angular) frequency, which gives the time rate of the phase 
evolution: 


UJ = 


lim <P(n) -^(0) _ 

n—> oo Ti 


( 4 ) 


III. NETWORKS OF COUPLED RULKOV NEURONS 

In the Rulkov model (JT j) -(|2 j) the variable x(n) plays the role of the membrane potential at 
discrete time t = nr, where r — 1 . Hence the difference x(n + 1 ) — x(n) can be interpreted 
as the time derivative of the potential. If the membrane capacitance is scaled to the unity, it 
amounts to the transmembrane current. Hence the effect of coupling is to inject a synaptic 
current in the equation for the x -variable ([I]). 

Let us denote by Xj the membrane potential of the ith neuron ( i = 1,2,... iV). The 
equations for a network of coupled neurons are 

N 

Of ’ X N 

Xi{n + 1) = -- - 2 + yi(n) + £ > A ijXj (n), (5) 

l+[Xi{n)] 

yi(n+l)=yi(n)-crxi(n)—p, ( 6 ) 

in which A t] is the adjacency matrix, which defines the type of synapse which connects the 
fth and jth neurons (respectively post-synaptic and pre-synaptic). 

If the spatial distance between them is taken into account, we can model an electrical 
synapse by a band-diagonal adjacency matrix which includes only near-neighbors of a given 
neuron. Chemical synapses are thus represented here by off-band-diagonal elements, since 
they include the effect of distant neurons. This is a simplified model, though, since it 
does not take into account the number of open channels in the post-synaptic neuron. For 
studies of synaptic-dependent phenomena, like plasticity, memory storage, learning, pattern 
recognition, etc. this simplified form of the coupling term may not be sufficient for numerical 
simulations. Moreover, we consider all interactions as bidirectional, in such a way that the 
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adjacency matrix is symmetric (Aj t = Aij), what is more likely to be the case in electrical 
than chemical synapses. However these simplifications are acceptable in a model as long as 
we are interested only in the effect of the network topology on the neuron dynamics. An 
example is the investigation of the DBS effects related to tremor associated with bursting 
synchronization in the thalamus 50]. 

The only parameter in our model that is to characterize the intensity of the synaptic 
connections is the coupling strength e. It takes on positive values for excitatory synapses 
and negative values for inhibitory ones. Since £ is the same for all neurons in our model we 
cannot consider here the case in which part of the synapses are excitatory (~ 75%) and part 
inhibitory. A modification such as this, although simple to implement in principle, would 
make the model more difficult to compare with minimal models of phase dynamics. The 
only restriction on the values of £ is that the coupling term itself cannot be too large so 
as to drive the neuron dynamics off a bursting state. This parameter must be adjusted by 
trial-and-error. 

Since in neuronal assemblies the neurons are likely to exhibit some diversity, we choose 
randomly the values of a within the interval [4.1,4.3] according to a specified probability 
distribution function (PDF) g(a) 51]. In spite of this we keep the other parameters with the 
same values (namely a = ft = 0.001) since this do not change our results in an appreciable 
way. We can compute the bursting phase and the corresponding angular frequencies 
a/*) using (J3]) and (|4j), respectively, for the coupled neurons in the same way as we did for 
isolated ones. 

When the a parameter is chosen in the interval [4.1,4.3] there is a linear relation between 
a and the mean burst frequency. For this reason when we choose the a parameter for each 
neuron according with a given PDF g{ot) we generate a PDF for the bursting frequency of 
uncoupled neurons denoted by g(Q) with the same shape as g{a). In this work we consider 
two types of distributions: (i) a waterbag (or uniform) distribution: 


„ for a < a < b 

~g w (a)={ a - b 

0 otherwise 


(7) 


where a = 4.1 and b = 4.3; and (ii) a truncated Cauchy distribution: 


g T (a) = 


7 


tan 


_i f b — cto 

7 


tan 


—i ( ^ 

7 


n -1 


1 + 


Q? — Q(q 

7 


-l 


( 8 ) 


8 











(a) 0 


T 


1 (b) Or 


-0.5 - - -0.5 

£3 -1 / p ^' -1 

>< -1-5 h -X-l.5 

-2 - 

41000 2 4(k)00 
(d) 2 




40500 


41000 


40500 


-2 (- 

41000 4&000 


H l Ml 111 1 


_i_ 

40500 


41000 



n 


n 


FIG. 2. Time evolution of (a) mean field and (c,e) fast variables of two selected neurons in a 
network of uncoupled Rulkov neurons. (b,d,f) stand for coupled neurons. 


where a = 4.1, b = 4.3, ato = 4.2 is the position of the distribution peak, and 7 = 0.1 is the 
half-width at half-maximum of the PDF 52], 

We shall defer the discussion of the nature and properties of the adjacency matrix to the 
following section. For the moment let us assume that a convenient form of has been 
given to the system of coupled Rulkov neurons flHJ-flHj)- One distinctive effect of synaptic 
coupling is bursting synchronization: two or more neurons begin their bursting activity at 
the same times (up to a given tolerance), regardless of whether or not the ensuing sequence 
of spikes coincides. A example is provided by Fig. [2] where we consider two randomly 
selected neurons in a network: in the uncoupled case the neurons burst at different and 
uncorrelated times [Fig. [21(c) and (e)], what makes the mean field 

1 N 

x (n) = jj ( 9 ) 

1=1 


to have small-amplitude fluctuations [Fig. [21(a)]. For large enough values of £ the coupled 
system of neurons present bursting synchronization as illustrated by the two selected neurons 
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[Fig. |2](d) and (f)]. The mean field in this case undergoes large-amplitude oscillations with 
the same frequency as of the bursting neurons [Fig. [2](b)]. 

Another quantitative characterization is provided by the order parameter magnitude and 


its time average given, respectively, by 

N 


33] 


R(n) = ^ 


i =1 


w>i(ri) 


R = -,Y. Rln >-- 


( 10 ) 


n =1 


where n' is chosen so as to yield stationary values of R. If the neurons burst exactly at the 
same time their phases coincide and hence the normalized sum in (ITU]) gives R — 1. On 
the other hand, if the neuron bursting is so uncorrelated that the times at which they burst 
are practically random the summation gives nearly zero and R ~ 0 for such an extreme 
case. In finite networks there is likely to be chance correlations, hence we consider R = 0.1 
as a threshold for the transition from a non-synchronized to a partially synchronized state. 
Intermediate cases (0.1 < R < 1) thus represent partially synchronized states. 

IV. CONNECTION ARCHITECTURES OF THE NEURONAL NETWORK 


Neurons are connected by axons and dendrites, so that we can regard those neurons as 
embedded in a three-dimensional lattice. However, due to the high connectivity of the neu¬ 
rons it is necessary to use complex networks to describe the properties of neuronal assemblies. 
These complex networks can be viewed in two basic levels: a microscopic, neuroanatomic 
level, and a macroscopic, functional level. Studies in the former level are limited to those 
few examples in which there is available data on the neuronal connectivity, as the worm C. 
Elegans [59]. 

In the second (macroscopic) level of description of neural networks, the use of non-invasive 
techniques as electroencephalography, magnetoencephalography and functional magnetic 
resonance imaging provides anatomical and functional connectivity patterns between dif¬ 


ferent brain areas 


60 


61]. This information provides a way to study the brain cortex, 


considering the latter as being divided into anatomic and functional areas, linked by axonal 
fibers. Scanned and coworkers have investigated the anatomical connectivity matrix of the 
visual cortex for the macaque monkey and the cat 62, [63]. 

The basic elements of a complex network are nodes and links. A link connects two or 
more nodes, and these connections can be unidirectional or bidirectional. In the language 
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of neuronal networks a node can be either a single neuron or a cortical area, depending on 
the level of characterization. If we consider a microscopic description the links are synapses 
(electrical or chemical), whereas they are axonal fibers in networks of cortical areas or simply 
they stand for the functional relationship between the corresponding regions in the cortex. 
A complex network can be characterized by various topological and metrical properties 


64]- For our purposes we will need only two of them: the average path length L and the 


average clustering coefficient C. The path length between two nodes a and b is the minimum 
number of links that must be traversed through the network to travel from a to b. The 
average path length L is the path length averaged over all pairs of nodes in the network. The 
average clustering coefficient C of a network is defined as the average fractions of pairs of 


neighbors of a given node which happen also to be neighbors of each other 


65]. For example, 


if a node a is connected to nodes b and c, and if b and c are themselves connected, then 
a, b and c form a triangular cluster. The value of C turns to be the fraction of triangular 
clusters that exist in the network with respect to the total possible number of such clusters. 

From the graph-theoretical perspective a complex network can be described by an ad¬ 
jacency matrix A t j , whose elements are equal to the unity if the nodes are connected and 
zero otherwise. In studies of cortical area networks, both in the anatomic and functional 
levels, each link is sometimes given a specific weight, which can be even time-varying in 
investigations of plasticity. Since we discard self-interactions, the diagonal elements of the 
adjacency matrix are zero. 


Global and random networks 


In a globally coupled network all nodes are connected with all other nodes. Since in this 
case the contribution of the coupling would increase with the number of nodes, the coupling 
term is usually divided by the number of nodes. In such a network each node can be regarded 
as being coupled to the mean held of other nodes. The adjacency matrix of global networks 
is constant: A* = 1 for a.. i,J = 1,2,...AT. There have been numerical investigations of 
such networks with respect to strategies of DBS [571.158]. The coupling term in ([5]) will be 
divided by the number of neurons, hence we shall redefine it as 


£ = 


1 

N 


( 11 ) 
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FIG. 3. (color online) Variation of the order parameter magnitude with the coupling strength for 
a globally coupled (a), and a randomly coupled (b) network of N = 1000 Rulkov neurons and 
probability p = 0.01. Black and red lines stand, respectively, for a waterbag and a truncated 
Cauchy distribution of values of the a parameter in the [4.1 : 4.3]. The green curves represent a 
polynomial fit given by Eq. (1161) . The insets zoom the behavior near the transition to bursting 
synchronization. 

The variation of the order parameter magnitude of the coupled network with the coupling 
strength above is depicted in Fig. [3](a), where we plot R as a function of £ when the 
values taken on by the parameter are drawn from a waterbag PDF (black line) and a 
truncated Cauchy distribution (red line). Since we used randomly chosen initial conditions 
for the neurons (x^, y^) we considered averages over 100 different realizations of the initial 
pattern. 

In both cases the behavior is qualitatively similar: for small coupling values there is no 
bursting synchronization and, after a critical coupling value £ c ~ 0.020 we begin to observe 
partial synchronization and, further on, complete phase synchronization. The approach to 
the latter differs with respect to the PDF g(a), being slightly faster for a truncated Cauchy 
distribution than for a waterbag PDF [see the inset of Fig. [31(a)]. 

The “physical” distance between two neurons does not play any role in globally coupled 
networks, since the contribution is counted in the same way for all nodes. This deficiency has 
been circumvented by considering a network in which the coupling strength decreases with 
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the lattice distance as a power-law, with exponent [70]. If q —* 0 the lattice distance does 
not matter and we recover the globally coupled case. As the exponent tends to infinity 
we approach a locally coupled network, where only the nearest-neighbor nodes have to be 
taken into account. 

Unlike global networks, random network present a typically small density of links. Ran¬ 
dom (Erdos-Renyi) networks are obtained by from A initially uncoupled nodes and build¬ 
ing Nk links between randomly chosen pairs with a uniform probability p (avoiding self¬ 
interactions) 66]. Since the total possible number of links is A(A — 1) one has 


N k 

,V(.V - 1)' 


( 12 ) 


The average path length of ER networks is typically very small since it scales as the logarithm 
of the network size ( L ran dom ~ In A”). Moreover the clustering coefficient of ER networks is 
likewise small since C ran dom = 1/A. 

Since the contribution of the coupling term is typically small we do not need here to rescale 
it as we did before. We show, in Fig. [3](b), the variation of the order parameter magnitude 
with the coupling strength e, showing a qualitatively similar picture but with a considerably 
lower critical coupling strength e c f« 0.002 for the transition from a non-synchronized to 
a partially synchronized bursting. If the probability p is too small there may not occur 
bursting synchronization at all, and we require a minimum number of links to observe such 
phenomenon, as shown in Fig. |3](b), where p = 0.01 and we have pN(N — 1) = 9990 links. 


B. Small world networks 

Random networks of the Erdos-Renyi type represent one limit case of a spectrum of 
networks of which the other end comprises regular networks, which are lattices of A nodes 
for which each node has links to its z nearest and next-to-the-nearest neighbors, hence there 
are local connections only. The average path length of one-dimensional regular lattices is 
Lreguiar ~ N and the clustering coefficient is given by C regu i ar = [3(2 — 2)]/[4(z — 1)] jfS?]]. 

In between those limiting cases we have small-world networks, for which we typically 
have L > L ran dom and C C ran dom ■ It is possible to obtain them from regular lattices 
essentially by two procedures: (i) Watts-Strogatz (WS) and (ii) Newman-Watts (NW) ones. 
WS networks are obtained from regular lattices by going through each of its links and, with 
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some specified probability p, randomly rewiring that link, moving one of its ends to a new 


node randomly chosen from the rest of the lattice 


68 ]. 


The coordination number of the network is still z on average, but the number of links per 
node may be greater or smaller than z due to the possible existence of nonlocal shortcuts, 
in addition to the local links. One possible disadvantage of this construction is that, if N 
is small enough, the rewiring process may create clusters disconnected to the rest of the 
network, for which nodes the path length would be obviously infinite. The latter problem 
can be circumvented in NW networks, that are constructed similarly to the WS ones, but 
the nonlocal shortcuts are added to the lattice with probability p, instead of being rewired 


69]. 


A small-world network of either WS or NW types is essentially described by the proba¬ 
bility of nonlocal shortcuts p. By computing both L and C we can get a range of p for which 
the small-world property is fulfilled by the network. The small-world property requires the 
network to exhibit a relatively small path length while retaining an appreciable degree of 
clustering, i.e. the following conditions L > L ranc ] or n and C C ran dom may hold for many 
different types of networks of the WS or NW types [68]. The adjacency matrix for such net¬ 
works is symmetric and band diagonal (with zero diagonal elements) and the non-diagonal 
parts are sparse, most of their elements being equal to zero. 

If the probability of shortcuts is zero we have the average path length and cluster co¬ 
efficient of a regular networks, denoted as L(0) and C'(O), respectively, for a network of 
N = 1000 neurons with z = 20 local connections. In Fig. [4](a) we plot the ratios C(p)/C( 0) 
and L(jp)/L( 0) as a function of the probability of randomly chosen shortcuts p in the NW 
procedure. The small-world property holds as long as the ratio C(p)/C( 0) is large and 
L(p)/L( 0) is small, what yields an interval [0.01, 0.1]. A further confirmation of this interval 
is provided by the merit figure <r, defined by 


K 


a= y 


where 


A = 


L(p) 


K, = 


Cjp) 

C T! 


(13) 


(14) 


-^rand '-''rand 

in such a way that the larger is a the better the small-world property holds for the network. 
We plot the merit figure as a function of p in Fig. Q](b), which indeed assumes larger values 
in the interval [0.01, 0.1]. The variation of the order parameter magnitude for a small-world 
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FIG. 4. (color online) (a) Dependence of the ratios C(p)/C( 0) (black line) and L(p)/L( 0) (red 
line) with the probability of shortcuts in a small-world network of the Newman-Watts type with 
N = 1000 and z = 20; (b) Dependence of the merit figure a with p; (c) Variation of the order 
parameter magnitude with the coupling strength for a small-world network of Rulkov neurons with 
probability p = 0.1. Black and red lines stand, respectively, for a waterbag and a truncated Cauchy 
distribution of values of the a parameter in the [4.1,4.3]. The green curve represents a polynomial 
fit given by Eq. (1161) . The insets zoom the behavior near the transition to bursting synchronization. 


network with p — 0.1 is depicted in Fig. [4](c) as a function of e, showing a transition to 
partial bursting synchronization for e c 0.001. 
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C. Scale-free networks 


In scale-free networks the connectivity ki (the number of links per node i) satisfies a 
power-law PDF 

P(k) ~ k-\ (7 > 1), (15) 

in such a way that highly connected neurons are connected, on the average, with highly 
connected ones, a property also found in many social and computer networks [64|, [74j]. Func¬ 
tional magnetic resonance imaging experiments have suggested that some brain activities 
can be assigned to scale-free networks, with a scaling exponent 7 between 2.0 and 2 . 2 , with 


a mean connectivity < k >~ 4 


71]. In fact, this scale-free property is consistent with the 


fact that the brain network increases its size by the addition of new neurons, and the latter 
attach preferentially to already well-connected neurons. 

In this paper we use the Barabasi-Albert (BA) coupling prescription through a sequence 
of steps, starting from an initial random (Erdos-Renyi) network of size N 0 = 23 nodes and 
23 random connections [3]. Every step we add a new node to the network which makes 
two connections: the first is determined by a uniform probability of connection among the 
vertices in the network and the second link is chosen such that the probability of connection 
decays according with the degree of connectivity of each node. Hence this second link is more 
likely to be attached to the most connected node in the network than to the less connected 
one. When the network reaches the desired size N we stop adding new nodes. 

We applied the BA procedure until the final network has N = 1000 nodes. We then 
obtained a numerical approximation for the PDF of connections per node P(k ), depicted 
in Fig. [5](a). The solid red line is a least squares fit giving a power-law dependence of 
the form (ITS]) with an exponent 7 = 2.9, confirming that the network is scale-free indeed. 
The variation of the order parameter magnitude with e is plotted in Fig. G3(b), illustrating 
the transition to bursting synchronization with e c = 0.004. The behavior of the average 
order parameter magnitude with e after the transition to bursting synchronization at e c 
can be fitted, when the frequency distribution is a truncated Cauchy PDF, by the following 
expression 


R 


1 - 


p \ 7~ 
C C 

£ 


( 16 ) 


where the exponents are different according to the type of network [Table H] . 
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FIG. 5. (color online) (a) Probability distribution function for the number of connections per 
site in a network obtained from BA procedure and N = 1000 nodes; (b) Variation of the order 
parameter magnitude with the coupling strength for a small-world network of Rulkov neurons with 
probability p = 0.1. Black and red lines stand, respectively, for a waterbag and a truncated Cauchy 
distribution of values of the a parameter in the [4.1 : 4.3]. The green curve represents a polynomial 
fit given by Eq. (1161) . The insets zoom the behavior near the transition to bursting synchronization. 


type 


r 

s 

Obs. 

global 

0.016/A 

4.5 

1.0 


Erdos-Renyi 

0.0017 

2.0 

1.0 

K = 10000 links 

small-world 

0.00075 

4.0 

2.0 

NW: p = 0.1, £ = 20 

scale-free 

0.004 

2.0 

0.7 

BA: No = 23, 7 = 2.9 


TABLE I. Critical coupling strength and fitting parameters for some types of networks with N = 
1000 nodes and a truncated Cauchy distribution. 

V. GENERALIZED KURAMOTO MODEL 

We can summarize the results of the previous Section by stating that, if the coupling 
strength is large enough in a network of coupled Rulkov neurons, there will be a transition 
to bursting synchronization when £ = e c . The critical values of the coupling strength are 
different, however, depending on the type of network (keeping the network size constant, see 
Table Q]) . This leads to the question of how the type of network might determine the value 
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of £ c . 


The answer to this question is elusive even for minimal models like the network of Rulkov 
maps considered in this work. However it can be formulated in a simpler model which has 
the virtue of being amenable to analytical methods like mean field theory, and which is the 
generalized Kuramoto model. Let us define a phase 0 t E [0, 2i r) for the All member of an 
assembly of N oscillators connected by a network of which we know the adjacency matrix 
Aj 

d6i 

dt 


N 


-Jj- = Ui + V ^2 A] sin (0j 

3 = 1 


ft), (i = l..... at) 


(17) 


where a is the coupling strength and oq are natural frequencies which we randomly choose 
from a PDF g(uj) which we require to be unimodal and symmetric: g(—u ;) = g(oj). The 
generalized Kuramoto model (GKM) has been used to study oscillations in cortical circuits 
75], as well as properties as axonal delay and synaptic plasticity 0 . 


A number of recent works has considered analytically the onset of synchronization in the 


GKM on complex networks 
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77]. Two basic approximations have been made: (i) the 


adjacency matrix is symmetric (A Jt = A t] ); (ii) the number of connections per node is large 
enough (/q 1). In this case the critical coupling is given by 

K r 


°cl 


Ar 


(18) 


where A r 


is the largest eigenvalue of the adjacency matrix and 

2 


K c = r \ 

*g( o) 


(19) 


is the critical coupling strength of the classical Kuramoto model, which corresponds to the 
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case of a globally coupled network 
A max = N — 1 and, at criticality, we have 


35]. Indeed, since A l3 = 1 for them it turns out that 


N 


dOi K c ^ . 

— = Ui H-> An sin 

dt. 1 n - y ■> 


6i). 


( 20 ) 


3 =1 


A mean field analysis based on the same assumptions as above furnishes the following 


estimate for the critical coupling strength in the GKM jr3] 

< k > 


c2 A r 


< k 2 >’ 


( 21 ) 
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where < k > and < k 2 > are the mean connectivity and variance, respectively, or the two 
first moments of the PDF P(k). The limitations of this formula are clear, though, since 
it was derived by assuming that hi 1, which is hardly the case for small networks. An 
example of the possible inadequacy of this formula is the case of a scale-free network, for 
which the variance 

/ oo />oo 

dkP{k)k 2 ~ J dkk 2 ~\ ( 22 ) 

which diverges if 7 < 3, yielding an infinite value for the critical coupling strength. 


VI. DISCUSSION 


In this Section we will compare results obtained from a network of coupled Rulkov net¬ 
works and the GKM. Similarities between them have been previously observed but no direct 


connection between those models has been made so far 


78 


80]. The starting point of our 


discussion is that the Kuramoto phase 9i(t) can be identified with the bursting phase of a 
Rulkov neuron If the time discretization step r is small enough it is immaterial if the 

model uses continuous time (like the GKM) or discrete time (like the Rulkov map). The 
similarity between these models can be proved on general grounds, when one is close to a 
globally phase synchronized state, the proof being sketched in an Appendix. 

The bursting frequencies, which are time rates of the corresponding phases, should be 
similar to both Rulkov and Kuramoto networks. Hence we adjust the PDF’s of the frequen¬ 
cies of the Rulkov model to comply with the symmetry requirements of the GKM. Let us 
first consider the waterbag distribution with b = —a: 


9w{u) — 


Ta for 

0 otherwise 


—a < ui < a 


(23) 


for which gw( 0) = l/2a. The parameter a for the GKM will be chosen so as to yield the 
same critical coupling strength which we numerically determined for the Rulkov network. 
For example, using the theoretical estimate (ITS]) we have at the critical point 


or, using CD 


_ 

°c 1 


\n 


7 r 5 , m(0)A 

r 


7r\ r 


(w) 2 < k > 

Gc2 ngwiti) < k 2 > 


Aa < k > 
7 t < k 2 > 


(24) 


(25) 
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network 

<k> 

< k 2 > 

^max 

global 

999 

999 2 

999 

Erdos-Renyi 

10 

109.30 

11.019 

small-world 

24.29 

594.83 

24.514 

scale-free 

3.954 

25.058 

6.33 


TABLE II. Properties of the network models used in this paper. 


It is possible to repeat the arguments for the case of a truncated Cauchy distribution 
with b = — a, a 0 = 0 , and where we have set co 0 = 0 and 7 = a in such a way that we have 

-1 


9t{u) = — 


Tia 


i + C 

a 


(26) 


with gr( 0) = 2 /( 7 ra). This yields theoretical estimates for the critical point 


a 


(T) 

cl 


a 


(T) 

c2 


2 a 

ng T (0)X max ^max 

2 < k > _ <k> 

7t<7t(0) < k 2 > a < k 2 > 


(27) 

(28) 


Let us consider a network of N — 1000 Rulkov neurons with the following connection 
topologies: global, random (ER), small-world (NW), and scale-free (BA), whose parameters 
are listed in Table QT| Using the numerically computed value of the critical coupling strength 
e c in the place of the a c in Eqs. (I24lh(l25]l and (I27|b(l28l) we can estimate the value of the 
parameter a of the PDF of natural frequencies using both the waterbag and the Truncated 
Cauchy models, respectively [cf. Tab. llllj . In the case of a random (ER) network, for 
example, in which e c = 0.020 and 0.016 for the waterbag and Cauchy PDFs, respectively, 
we found a = 0.016 and 0.017. Similar results hold for other network topologies. The values 
of a c 1 and a C 2 present only slight differences. 

After having chosen adequate values of the parameters appearing in the PDF of both 
Rulkov and Kuramoto models, we can compare the results for them with respect to the 
transition to bursting synchronization. I 11 Figs. El(a) to (d) we compute the order parameter 
magnitude for the Rulkov (black line: waterbag, red line: Cauchy) and Kuramoto (green 
line: waterbag, blue line: Cauchy) as a function of the coupling strength for the network 
topologies considered in this paper. The evolution of R is nearly the same for both models, 
with small differences. For example the value of R in the GKM reaches its maximum slightly 
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network 

waterbag 

truncated Cauchy 

global 

0.016 

0.017 

Erdos-Renyi 

0.017 

0.019 

small-world 

0.017 

0.018 

scale-free 

0.025 

0.025 


TABLE III. Estimates for the value of parameter a for two different PDFs and some types of 
networks. 

before (i.e. with smaller coupling strengths) the Rulkov network does. This precedence of 
GKM in comparison with Rulkov is observed for all connection topologies we have considered 
in this work. 

VII. CONCLUSIONS 

In this work we propose minimal models for the description of bursting neurons using 
some of the most used connection topologies for neuronal networks which can be used in nu¬ 
merical experiments. One of the simplest models that describe bursting is a two-dimensional 
discrete-time mapping proposed by Rulkov, and we considered networks of coupled Rulkov 
neurons using global, random (Erdos-Renyi), small-world (Newman-Watts), and scale-free 
(Barabasi-Albert) topologies. We can define a geometrical phase for coupled Rulkov neu¬ 
rons, such that bursting synchronization is equivalent to some form of phase synchronization. 
Using the time rate of the bursting phase we can define a bursting frequency and thus in¬ 
vestigate frequency synchronization in such networks. 

Using a complex order parameter we quantify the state of phase synchronization of the 
system. As the coupling strength is increased, there is a transition from a non-synchronized 
to a partially phase-synchronized state, which evolves to a completely phase-synchronized 
one for larger values of the coupling strength. We computed the corresponding critical values 
of the coupling strength for the connection topologies considered in this work. 

The existence of a phase reduction for the coupled Rulkov neuron system suggests the 
usefulness of a Generalized Kuramoto model which, in the present context, can be considered 
a minimal model for bursting neurons. We presented in this paper some quantitative results 
that justify this claim. In order to compare both models, however, it is necessary to adjust 
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8 8 


FIG. 6. (color online) Order parameter magnitude vs. coupling strength for a network of Rulkov 
neurons (black: waterbag, red: Cauchy) and Kuramoto oscillators (green: waterbag, blue: Cauchy) 
with N = 1000 nodes, for (a) global coupling, (b) Erdos-Renyi; (c) small-world, and (d) scale-free. 


the corresponding probability distribution functions for the natural frequencies. We did so 
by using theoretical estimates for the critical coupling strength for the Generalized Kuramoto 
model in different connection topologies. After choosing the PDF parameters in a suitable 
way the behavior of the order parameter is very similar for both the Rulkov and Kuramoto 
models. 


22 






































Appendix A: Phase reduction near global phase synchronization 


In this appendix we show that, in the scenario near a phase synchronized regime, one 
can perform a phase reduction from a network of coupled oscillators to a generalized Ku- 
rarnoto model, what explains why the properties of a coupled Rulkov neuronal network are 
similar to the Kuramoto model. We stress that this similarity does not necessarily apply to 
the situation near a frequency synchronized regime, since the dynamical requirements are 
different. Our presentation follows closely that in Ref. 79J and the general treatment given 


in Ref. 


33]. 


Let a general .D-dimensional flow be given by 

dx 


dt 


— Fix') 


(Al) 


where x is a D-dimensional vector in the system phase space and F a vector held. We 
assume that there is a stable period-T orbit 


xo(t) = x 0 (f + T). (A2) 

The “slow” dynamics along this periodic orbit can be described by a phase cp(x) which time 
evolution is given by 

-^(x) = V x ^F(x) = 1. (A3) 

Now let us consider a network of coupled oscillators with a slight mismatch in their 
parameters described by 

^Xj = F(xj) + fife) + ey^a < jV(x i ,Xj), (i = 1, 2,... N), (A4) 

3 

where is different for each oscillator and stands for the vector held part containing slightly 
mismatched parameters, e is the coupling strength, a t] the adjacency matrix, and V(x, : , x^) 
is a coupling function. 

From (I A 31) . the slow phase ipi of the coupled oscillators is implicitely defined by the 
function 

Z {ipi) = V x ^(xo(^i)), (A5) 

where xo is a period-T stable orbit. The time evolution of the phase is then governed by 

^x = 1 + ey^QijZ(y> i )V(x i ,Xj) + Z(<£>j)fj(xj), (A6) 

3 
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for i — 1, 2,... N. 

On introducing an auxiliary phase = fi — 1 and using a time average over a period T, 
by keeping ?/;,■ fixed we can write, in a first-order approximation, an equation governing the 
time evolution of vjjf 


d 


—^ = iDj + e ayT^i, ipj), 


(AT) 


where 


w; = 


T 

1 


Z(t + 


= — / z(t + V’t)v(x 0 (t + v , i) ; x o(^ + i>j))dt, 


(AS) 

(A9) 


play the roles of frequencies and coupling functions, respectively, for the auxiliary phases 
If we are close to a global phase synchronized state, for which ipi ps (pj for any pairs of 
oscillators (i,j), we have an approximate form for the coupling function T. Introducing the 
time-dependent variable ( = t + ipj and supposing that ^ T there results 


X fT+'Pi 

r tyiM = f J z(c + '0i-'0i)v(x o (C + '0i-'0j), x o(C + '0j))^- (Aio) 

By expanding the coupling function V in a power series and assuming that Z is nearly 
constant over the periodic orbit x 0 , we get 


rOi, i>j) ~ a + b(ipi - ijjj), 


where 

« - ^ J ZV(xo(C),xo(C))<iCi 

Z(V x V(x,y)) x _ y=xo(0 ^J(, 

Notice that this approximation holds whenever 

V x V(x,y)^0 


for any x and y belonging to the periodic orbit xo- 
Substituting (IA11D into (IA7D yields 



Wi + Sid + eb^2aij(ipi —■ -0j), 

j 


(All) 


(A12) 

(M3) 


(A 14) 


(A 15) 
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where s t = JA a ij is the intensity of the ith node of the network. As long as we deal with 
networks for which the intensities s* present only a small variation over the network, we can 
take Si as practically constant, i.e., independent of i. 

Moreover, since we are by hypothesis near a global phase synchronized state, the phase 
difference ipi — t/ij is small enough to justify the replacement ipi — ~ sin('0 i — ifj), in 

such a way that the equation governing the time evolution of the auxiliary phases (near a 
phase-synchronized situation) is a generalized Kuramoto model 

ipi^Wj + e a i:j sin- ipj), (A16) 

j 

where Wi = Wi + asi and £ = eb. 
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